
# Run analysis pilot ------------------------------------------------------

fire_teach_pilot <- ols_main(
	outcome = "fire_teach",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)


fire_teach_pilot_pvals <- get_RI_pvals(
	outcome = "fire_teach",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment_pilot,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

pta_pilot <- ols_main(
	outcome = "pta",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)


pta_pilot_pvals <- get_RI_pvals(
	outcome = "pta",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment_pilot,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

bring_up_pilot <- ols_main(
	outcome = "bring_up",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)


bring_up_pilot_pvals <- get_RI_pvals(
	outcome = "bring_up",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment_pilot,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

assemble_pilot <- ols_main(
	outcome = "assemble",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)


assemble_pilot_pvals <- get_RI_pvals(
	outcome = "assemble",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment_pilot,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

absenteeism_action_pilot <- ols_main(
	outcome = "absenteeism_action",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

absenteeism_action_pilot_pvals <- get_RI_pvals(
	outcome = "absenteeism_action",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(el_pilot, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment_pilot,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")


# Run analysis main experiment midline ----------------------------------------------------

lc1_ml <- ols_main(
	outcome = "lc1",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

lc1_ml_pvals <- get_RI_pvals(
	outcome = "lc1",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

pta_ml <- ols_main(
	outcome = "pta",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

pta_ml_pvals <- get_RI_pvals(
	outcome = "pta",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

bring_up_ml <- ols_main(
	outcome = "bring_up",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

bring_up_ml_pvals <- get_RI_pvals(
	outcome = "bring_up",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")

assemble_ml <- ols_main(
	outcome = "assemble",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

assemble_ml_pvals <- get_RI_pvals(
	outcome = "assemble",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")


absenteeism_action_ml <- ols_main(
	outcome = "absenteeism_action",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE)

absenteeism_action_ml_pvals <- get_RI_pvals(
	outcome = "absenteeism_action",
	treatment = "absenteeism",
	resample_FE = TRUE,
	block_FE = TRUE,
	audience_size = TRUE,
	cluster_SE = TRUE,
	covariates = NULL,
	the_data = subset(ml, compliance == 1),
	dosage = FALSE,
	dosage_indicator = FALSE,
	assignment_data = treatment_assignment,
	extract_function = coef,
	analysis_function = ols_main,
	sims = sims,
	lwr_upr_two = "upr")



#Running Bayesian meta-analysis to combine pilot and main study results -------------------------

##Perform Bayesian meta-analysis with uninformative priors
#create handy function
meta_analysis <- function(
	pilot_b,
	pilot_se,
	main_b,
	main_se,
	...
){
	pilot_wt = 1/pilot_se^2
	main_wt = 1/main_se^2
	sum_wt = pilot_wt + main_wt
	
	pilot_estwt = pilot_wt * pilot_b
	main_estwt = main_wt * main_b
	sum_estwt = pilot_estwt + main_estwt
	
	overall_b = sum_estwt / sum_wt
	overall_se = 1/sqrt(sum_wt)
	overall_t = abs(overall_b / overall_se) #computes absolute value
	overall_p = 1 - pnorm(overall_t)

	output = c(overall_b,overall_se,overall_t,overall_p)
	names(output) <- c("Estimate","Std. Error","t value", "Pr(>|t|)")
	return(output)
	
}


#run meta-analyses - requires results from pilot_field_experiment_analysis.R and conative_attitudes_compliers.R

fire_teach_meta <- meta_analysis(fire_teach_pilot$fit_summary["absenteeism","Estimate"],fire_teach_pilot$fit_summary["absenteeism","Std. Error"],
																 lc1_ml$fit_summary["absenteeism","Estimate"],lc1_ml$fit_summary["absenteeism","Std. Error"],)

bring_up_meta <- meta_analysis(bring_up_pilot$fit_summary["absenteeism","Estimate"],bring_up_pilot$fit_summary["absenteeism","Std. Error"],
															 bring_up_ml$fit_summary["absenteeism","Estimate"],bring_up_ml$fit_summary["absenteeism","Std. Error"],)

pta_meta <- meta_analysis(pta_pilot$fit_summary["absenteeism","Estimate"],pta_pilot$fit_summary["absenteeism","Std. Error"],
															 pta_ml$fit_summary["absenteeism","Estimate"],pta_ml$fit_summary["absenteeism","Std. Error"],)

assemble_meta <- meta_analysis(assemble_pilot$fit_summary["absenteeism","Estimate"],assemble_pilot$fit_summary["absenteeism","Std. Error"],
															 assemble_ml$fit_summary["absenteeism","Estimate"],assemble_ml$fit_summary["absenteeism","Std. Error"],)

absenteeism_action_meta <- meta_analysis(absenteeism_action_pilot$fit_summary["absenteeism","Estimate"],absenteeism_action_pilot$fit_summary["absenteeism","Std. Error"],
															 absenteeism_action_ml$fit_summary["absenteeism","Estimate"],absenteeism_action_ml$fit_summary["absenteeism","Std. Error"],)


##Rather clunky way to get the tables to output properly in Stargazer
fire_teach_meta_coeftest <- fire_teach_pilot #duplicate coeftest object as a container for the meta-analysis results
fire_teach_meta_coeftest$fit$coefficients[] <- NA #zero out the statistics
fire_teach_meta_coeftest$fit_summary[,] <- NA #zero out the statistics
fire_teach_meta_coeftest$desc_stats[] <- NA #zero out the statistics
fire_teach_meta_coeftest #check that everything is NA
fire_teach_meta_coeftest$fit$coefficients["absenteeism"] <- fire_teach_meta["Estimate"] #put the b coefficient into the container
fire_teach_meta_coeftest$fit_summary["absenteeism",] <- fire_teach_meta #put the rest of the statistics in.  


bring_up_meta_coeftest <- bring_up_pilot #duplicate coeftest object as a container for the meta-analysis results
bring_up_meta_coeftest$fit$coefficients[] <- NA #zero out the statistics
bring_up_meta_coeftest$fit_summary[,] <- NA #zero out the statistics
bring_up_meta_coeftest$desc_stats[] <- NA #zero out the statistics
bring_up_meta_coeftest #check that everything is NA
bring_up_meta_coeftest$fit$coefficients["absenteeism"] <- bring_up_meta["Estimate"] #put the b coefficient into the container
bring_up_meta_coeftest$fit_summary["absenteeism",] <- bring_up_meta #put the rest of the statistics in.  

pta_meta_coeftest <- pta_pilot #duplicate coeftest object as a container for the meta-analysis results
pta_meta_coeftest$fit$coefficients[] <- NA #zero out the statistics
pta_meta_coeftest$fit_summary[,] <- NA #zero out the statistics
pta_meta_coeftest$desc_stats[] <- NA #zero out the statistics
pta_meta_coeftest #check that everything is NA
pta_meta_coeftest$fit$coefficients["absenteeism"] <- pta_meta["Estimate"] #put the b coefficient into the container
pta_meta_coeftest$fit_summary["absenteeism",] <- pta_meta #put the rest of the statistics in.  

assemble_meta_coeftest <- assemble_pilot #duplicate coeftest object as a container for the meta-analysis results
assemble_meta_coeftest$fit$coefficients[] <- NA #zero out the statistics
assemble_meta_coeftest$fit_summary[,] <- NA #zero out the statistics
assemble_meta_coeftest$desc_stats[] <- NA #zero out the statistics
assemble_meta_coeftest #check that everything is NA
assemble_meta_coeftest$fit$coefficients["absenteeism"] <- assemble_meta["Estimate"] #put the b coefficient into the container
assemble_meta_coeftest$fit_summary["absenteeism",] <- assemble_meta #put the rest of the statistics in.  

absenteeism_action_meta_coeftest <- absenteeism_action_pilot #duplicate coeftest object as a container for the meta-analysis results
absenteeism_action_meta_coeftest$fit$coefficients[] <- NA #zero out the statistics
absenteeism_action_meta_coeftest$fit_summary[,] <- NA #zero out the statistics
absenteeism_action_meta_coeftest$desc_stats[] <- NA #zero out the statistics
absenteeism_action_meta_coeftest #check that everything is NA
absenteeism_action_meta_coeftest$fit$coefficients["absenteeism"] <- absenteeism_action_meta["Estimate"] #put the b coefficient into the container
absenteeism_action_meta_coeftest$fit_summary["absenteeism",] <- absenteeism_action_meta #put the rest of the statistics in.  

#produce table  --------------------------------------------------
#setup

pval_lines <- c(
	"$p$-values",
	round(fire_teach_pilot_pvals$ri_pvals["absenteeism"], 2),
	round(lc1_ml_pvals$ri_pvals["absenteeism"], 2),
	round(fire_teach_meta_coeftest$fit_summary["absenteeism","Pr(>|t|)"], 2),
	round(bring_up_pilot_pvals$ri_pvals["absenteeism"],2),
	round(bring_up_ml_pvals$ri_pvals["absenteeism"],2),
	round(bring_up_meta_coeftest$fit_summary["absenteeism","Pr(>|t|)"],2),
	round(pta_pilot_pvals$ri_pvals["absenteeism"],2),
	round(pta_ml_pvals$ri_pvals["absenteeism"],2),
	round(pta_meta_coeftest$fit_summary["absenteeism","Pr(>|t|)"],2),
	round(assemble_pilot_pvals$ri_pvals["absenteeism"],2),
	round(assemble_ml_pvals$ri_pvals["absenteeism"],2),
	round(assemble_meta_coeftest$fit_summary["absenteeism","Pr(>|t|)"],2)
)

hypothesis_lines <- c(
	"Hypothesis",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr",
	"upr"
)

#make tables
sink("03_tables/ABS_conative_meta_analysis.tex")
stargazer(
	... = list(
		fire_teach_pilot$fit, 
		lc1_ml$fit,
		fire_teach_meta_coeftest$fit,
		bring_up_pilot$fit,
		bring_up_ml$fit,
		bring_up_meta_coeftest$fit,
		pta_pilot$fit,
		pta_ml$fit,
		pta_meta_coeftest$fit,
		assemble_pilot$fit,
		assemble_ml$fit,
		assemble_meta_coeftest$fit
	),
	type = "latex",
	p = list(
		fire_teach_pilot_pvals$ri_pvals, 
		lc1_ml_pvals$ri_pvals,
		fire_teach_meta_coeftest$fit_summary[,"Pr(>|t|)"],
		bring_up_pilot_pvals$ri_pvals,
		bring_up_ml_pvals$ri_pvals,
		bring_up_meta_coeftest$fit_summary[,"Pr(>|t|)"],
		pta_pilot_pvals$ri_pvals,
		pta_ml_pvals$ri_pvals,
		pta_meta_coeftest$fit_summary[,"Pr(>|t|)"],
		assemble_pilot_pvals$ri_pvals,
		assemble_ml_pvals$ri_pvals,
		assemble_meta_coeftest$fit_summary[,"Pr(>|t|)"]	),
	se = list(
		fire_teach_pilot$fit_summary[,"Std. Error"], 
		lc1_ml$fit_summary[,"Std. Error"],
		fire_teach_meta_coeftest$fit_summary[,"Std. Error"],
		bring_up_pilot$fit_summary[,"Std. Error"],
		bring_up_ml$fit_summary[,"Std. Error"],
		bring_up_meta_coeftest$fit_summary[,"Std. Error"],
		pta_pilot$fit_summary[,"Std. Error"],
		pta_ml$fit_summary[,"Std. Error"],
		pta_meta_coeftest$fit_summary[,"Std. Error"],
		assemble_pilot$fit_summary[,"Std. Error"],
		assemble_ml$fit_summary[,"Std. Error"],
		assemble_meta_coeftest$fit_summary[,"Std. Error"]
	),
	keep = "absenteeism",
	omit.stat = c("adj.rsq","n","rsq","f","ser"), #dropping the n and adj rsquared
	column.separate = c(3,3,3,3),
	column.labels = c("Ask LC1/Headmaster","Tell village","Use PTA","Assemble group"),
	table.layout = "=cd#-t-as=n",
	dep.var.labels = c("Pilot","Main","Meta",
										 "Pilot","Main","Meta",
										 "Pilot","Main","Meta",
										 "Pilot","Main","Meta"
	),
	dep.var.labels.include = TRUE,
	no.space = T,
	omit = "block_id",
	add.lines = list(
		pval_lines,
		hypothesis_lines,
		c("Observations",
			"376","1,041","1,417",
			"376","1,041","1,417",
			"376","1,041","1,417",
			"376","1,041","1,417")), #manually entering the n
	notes.label = "",
	float = FALSE,
	multicolumn = F
	# style = "qje"
)
sink()



pval_lines <- c(
	"$p$-values",
	round(absenteeism_action_pilot_pvals$ri_pvals["absenteeism"],3), 
	round(absenteeism_action_ml_pvals$ri_pvals["absenteeism"],3),
	round(absenteeism_action_meta_coeftest$fit_summary["absenteeism","Pr(>|t|)"],3)
)

hypothesis_lines <- c(
	"Hypothesis",
	"upr",
	"upr",
	"upr"
)


#make tables
sink("03_tables/ABS_conative_meta_analysis_index.tex")
stargazer(
	... = list(
		absenteeism_action_pilot$fit,
		absenteeism_action_ml$fit,
		absenteeism_action_meta_coeftest$fit
	),
	type = "latex",
	p = list(
		absenteeism_action_pilot_pvals$ri_pvals,
		absenteeism_action_ml_pvals$ri_pvals,
		absenteeism_action_meta_coeftest$fit_summary[,"Pr(>|t|)"]
	),
	se = list(
		absenteeism_action_pilot$fit_summary[,"Std. Error"],
		absenteeism_action_ml$fit_summary[,"Std. Error"],
		absenteeism_action_meta_coeftest$fit_summary[,"Std. Error"]
	),
	keep = "absenteeism",
	omit.stat = c("adj.rsq","n","rsq","f","ser"), #dropping the n and adj rsquare
	column.separate = c(3),
	column.labels = c("Conative attitudes index"),
	table.layout = "=cd#-t-as=n",
	dep.var.labels = c("Pilot","Main","Meta"
	),
	dep.var.labels.include = TRUE,
	no.space = T,
	omit = "block_id",
	add.lines = list(
		pval_lines,
		hypothesis_lines,
		c("Observations",
			"376","1,041","1,417")), #manually entering the n
	notes.label = "",
	float = FALSE,multicolumn = F
	# style = "qje"
)
sink()
